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Motivated by recent experiments, where a voltage biased Josephson junction is placed in series 
with a resonator, the classical dynamics of the circuit is studied in various domains of parameter 
space. This problem can be mapped onto the dissipative motion of a single degree of freedom in 
a nonlinear time-dependent potential, where in contrast to conventional settings the nonlinearity 
appears in the driving while the static potential is purely harmonic. For long times the system ap¬ 
proaches steady states which are analyzed in the underdamped regime over the full range of driving 
parameters including the fundamental resonance as well as higher and sub-harmonics. Observables 
such as the dc-Josephson current and the radiated microwave power give direct information about 
the underlying dynamics covering phenomena as bifurcations, irregular motion, up- and down con¬ 
version. Due to their tunability, present and future set-ups provide versatile platforms to explore 
the changeover from linear response to strongly nonlinear behavior in driven dissipative systems 
under well defined conditions. 
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I. INTRODUCTION 

The nonlinear properties of Josephson junctions (JJs) 
have made such devices a key circuit element for clas¬ 
sical and quantum electronics. Accordingly, there has 
been a long tradition of studying non-linear phenomena 
in driven superconducting circuits, starting as early as 
the 1960s with the discovery of Shapiro steps^. While 
Shapiro-steps have remained a tool in exploring new 
directions in Josephson physics, for instance in atomic 
point contacts^^®, other nonlinear phenomena, like syn¬ 
chronization, have been investigated in arrays of JJs: 
as a test-bed for generic theory models to capture syn¬ 
chronization phenomena^’®, but as importantly with the 
prospect of applications as sources of more intense coher¬ 
ent radiation®, cf. also new developments using intrinsic 
arrays^®^^®. 

More recently, the nonlinearity of the JJ was ex¬ 
ploited as the crucial factor in enabling the high sensi¬ 
tivity of Josephson bifurcation amplifiers, achieving sub¬ 
stantial improvements towards reaching quantum-limited 
measurement processes^®^^®. Most of the features of 
Josephson bifurcation amplifiers, in fact, only rely on 
(and can consequently be described by) any type of 
nonlinearitye.g.. Duffing-type models, so that only 
recently the full nonlinear potential of the JJ has become 
of interest in this field^®. 

A recent addition to the field of driven nonlinear 
JJs®^“®® are experiments on a dc-biased JJ connected 
to a resonator®^^®®. In this sort of setup, charge trans¬ 
fer through the JJ leads to excitations in the resonator, 
and therefore allows to convert energy carried by charge 
quanta into quantum microwave photons. In these de¬ 
vices measurement of both the Josephson current and the 
emitted microwave radiation is possible, a distinct advan¬ 
tage in comparison to other recently proposed transport 


setups®®’®^, which show similar nonlinear features like 
bifurcations, period multiplication and up- and down- 
conversion®®^®^. Such nonlinear effects will, in fact, dom¬ 
inate the system’s dynamics and therewith the charac¬ 
teristics of the Josephson current and the emitted mi¬ 
crowave radiation for driving beyond a linear regime of 
weak Josephson coupling. 

While quantum properties of the system, in particu¬ 
lar of the emitted microwave radiation have been inves¬ 
tigated widely, both theoretically®®^"'® and within new 
experiments^®, due to the complexity and richness of non¬ 
linear effects a deeper understanding of the purely clas¬ 
sical dynamics of the system is instructive, but also of 
relevance for current experimental activities®®, as we will 
discuss below. This is especially the case, as the nonlin¬ 
earity enters the system in a peculiar way. It does not 
stem from a nonlinear (static) potential, but rather from 
the manner of coupling the drive (JJ) to the resonator. 
The Josephson phase in this setup is thus not fixed by 
the external voltage, but appears as a dynamical degree 
of freedom manifest in a time-dependent effective poten¬ 
tial determining the phase dynamics. We note in passing 
that related nonlinear phenomena have recently been ex¬ 
plored for Josephson phase slip devices^^. 

In this paper we present analytical and numerical in¬ 
vestigations in the regime, where the system’s dynamics 
is described by classical Josephson equations. While the 
features found are to an extent common to a wide class of 
nonlinear classical systems, which specific effects in what 
distinct manner are realized and how they are observed in 
this new type of nonlinear system is an intriguing open 
question. To tackle it, we will first present the system 
under study in more detail and introduce the analytical 
methods used in Sec. II. The following sections cover the 
fundamental resonance (Sec. Ill) and higher order reso¬ 
nances (Sec. IV). The influence of a thermal environment 


at finite temperatures is investigated in Sec. V, before we 
conclude in Sec. VI. 


2 


II. CIRCUIT DYNAMICS 


We consider a circuit (see Fig. 1), where a Josephson 
junction (JJ) is placed in series with a resonator with only 
a single mode of frequency ujq = Ij'jLC being relevant. 
The total impedance Zt{uj) seen by the tunneling Cooper 
pairs consists of the combination of the capacitance Cj of 
the JJ and the parallel LC resonance with finite Q factor. 
Since experimentally Cj C, one has ~ Z(uj) 

with 

= /n ■ (1) 

G — luq) + ujijJqIQ 

Based on Kirchhoff rules, equations of motion for the 
circuit, biased by a dc-voltage V, are found, 


(p+ujg(p+u;o^ + (Ej/(/>QC) sin(</)) = 0 , </) = ip+ujjt , (2) 

expressed in terms of the resonator’s phase variable 
(p = —{2efh) f dtVLc(i), and the Josephson phase (j) = 
l2e/h) J dtVj{t) with V = Vlc + Here, wj = 

2eVITi denotes the driving frequency, (j)o = ?i/2e the re¬ 
duced flux quantum, and the time derivative • = d/dt. 

This set of equations can be cast into an equation of 
motion for a fictitious particle with an effective mass, 
m = iPqC, moving in a harmonic potential and coupled 
to an external time-periodic, position dependent force, 
i.e.. 


Wo 


rrvp J- mujQip J- m — ip J- Ej sin ((/3 J- ujt) = ^{t) . (3) 


To incorporate finite temperature effects, as discussed 
below, thermal noise f at temperature T is added.It 
is related to the resonator damping via the fluctuation- 
dissipation theorem as p = 2m{oJo/Q)kBT6{t — 

t') and ( 0/3 = 0 . 

To explore the dynamics of the above Langevin equa¬ 
tion, it is convenient to work with dimensionless units, 
where times are scaled with ujq and energies with mujQ. 
This then leads to 


ifi + + X sm{nt + (f) = ^{t), (4) 

where the dimensionless friction coefficient 7 is related 
to the Q-factor of the resonator via 7 = 1/Q and 
we further introduced the dimensionless driving ampli¬ 
tude A = Ej/mcog and driving frequency fl = ujjjujo. 
Throughout the first part of the paper we concentrate on 
the limit T = 0 (^ = ^jmujQ = 0) and discuss the impact 
of thermal noise later. 

This form of the equation of motion, (4) , is the starting 
point for studying the dynamics of the voltage-biased cir¬ 
cuit (cf. Fig. 1) in the rest of this paper. In particular, we 
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Figure 1. Circuit diagram of a voltage-biased Josephson junc¬ 
tion in series with a resonator with relevant mode frequency 
cuq described by an effective impedance Z(uj) as specified in 
(1). Using a SQUID-geometry^ for the JJ allows for tuning 
the effective Josephson energy Ej(^x) by an external mag¬ 
netic flux 4>a;. 


are interested in the long-time limit, where the balance 
between the dissipative and the driving part of Eq. (4) has 
pushed the system into time-periodic steady-state orbits. 
Considering the energy balance of the resonator obtained 
from Eq.(4), 

It * t) = 

— Fhiss T , 

we easily identify the power dissipated from the res¬ 
onator, Pdiss) and the power injected into the resonator 
via the driving, Fjj_>.ho = Ij^lc, which will both be 
considered in detail below. The dissipated power to¬ 
gether with the dimensionless Josephson current Ij = 
A sin((/j J- Ht) constitute the main observables, which can 
be accessed experimentally, either averaged over many 
oscillation periods, or time- or frequency resolved. Note 
that the physical current results from the dimensionless 
current by multiplying (po/E- 

Due to the nonlinearity present in (4) the structure of 
steady-state orbits will depend sensitively on the damp¬ 
ing and amplitude and frequency of the driving, giving 
rise to the full wealth of nonlinear phenomena such as 
bifurcations, up- and down-conversion etc. What dis¬ 
tinguishes the situation under consideration here from 
most other driven nonlinear systems is that the nonlin¬ 
earity appears not in form of a static potential energy but 
rather as part of the driving force. In fact, it turns out 
that the effective time-dependent potential giving rise to 
(4), namely. 


can sometimes be illuminating to achieve a better under¬ 
standing of the fictitious particle dynamics. 

Of course, for nonlinear, time-dependent problems 
such as (4) explicit solutions can, in general, only be ob¬ 
tained numerically. However, analytical progress, at least 
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driving strength 


Figure 2. (color online) Tuning the dimensionless Josephson coupling A the steady-state dynamics at the fundamental resonance, 
n = Wj/cuo = 1, accesses different ranges discussed in Secs. IIIA, IIIB, IIIC, which are characterized by typical phase space 
portraits. Within the domain III B, the dotted vertical line separates the range of regular dynamics from the range A 2 < A, where 
the transition towards irregular dynamics occurs. The solid line illustrates the corresponding behavior of the dimensionless 
dc-current through the JJ with the shadow indicating its standard deviation when starting the steady-state dynamics with 
thermal initial conditions. The specific data shown for 7 = 0.01 exemplify the generic behavior in the underdamped regime. 
In this regime, the parameters A 1 / 7 , A* and A 2 are not very sensitive to varying the friction strength. 


for the steady state, can be made in limiting domains of 
parameter space. Physically, one expects a steady-state 
solution oscillating with the frequency of the drive, as 
well as higher and possibly subharmonics. Putting such 
an ansatz into (4), the nonlinear sin-term will again pro¬ 
duce higher harmonics with the coefficients of the various 
frequency components appearing in Bessel-functions. 

Formally, one can write a generic ansatz, 

keMn 

where ipk = ^-k and 9k — —9-k so that solutions 
are real-valued. The sum runs over a suitable set of 
rational numbers Mq with the prime indicating that 
fc = 0 is excluded. For example, below we will study 
the fundamental resonance at = 1 implying Mq = Z. 
Other situations are driving at O = n (n integer) with 
Mn = {u/n\v e Z} and at O = 1/n (n integer) with 
Mq = Z. Of course, in a perturbative treatment only a 
finite number of coefficients is taken into account. 

Now, inserting this ansatz into (4) yields a nonlinear 
equation for the steady-state Fourier coefficients, i.e., 

( 8 ) 

+ A re*^‘e*‘^«T+(t)-e-*^‘e-*^«F_(t)l =0, 

2i ^ 


with A/j = — 1 and F±{t) = H F^{ipk,t). These 

k>0 

latter functions 

00 

I —— 00 

contain Bessel functions Ji{-) of the first kind and of in¬ 
teger order. Eq. ( 8 ) serves as a starting point for per¬ 
turbative treatments in various ranges of the dynamics 
in order to gain a deeper understanding of the numerical 
findings based on (4). 

III. FUNDAMENTAL RESONANCE 

We will start in this section with the dynamics near 
and at the fundamental resonance where the driving fre¬ 
quency matches the Josephson frequency so that 0 = 1. 
According to (4), in absence of noise only two dimen¬ 
sionless parameters are left which determine the nature 
of steady-state orbits, namely, the friction 7 and the driv¬ 
ing amplitude A. Current experimental realizations are 
operated in the underdamped regime with a fixed 7 <C 1 
(Q-factors vary from 10 to about 1000 ) and varying driv¬ 
ing strengths. In the classical domain that we consider 
here, this leads from simple linear to strongly nonlinear 
phase space patterns, the structure of which is reflected 
in specific observables such as charge current and photon 
flux. 
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Before we study the details, let us give a brief qualita¬ 
tive account of the different dynamical ranges when A is 
varied (cf. Fig. 2): 

(i) In the regime of weak driving (Sec. Ill A), the dynam¬ 
ics changes from linear to nonlinear towards a threshold 
Ai, where a hrst bifurcation occurs. Perturbative treat¬ 
ments capture this transition fairly accurate. While the 
oscillation amplitude already reflects nonlinearities of the 
system, phase-space orbits in this regime are basically 
ellipses with only small deviations from the harmonic 
limit. Physically, the initial quadratic rise of the dc- 
current through the JJ flattens with increasing A until 
it saturates at A = Ai. 

(ii) For an intermediate range Ai < A < A 2 (Sec. IIIB) 
the dynamics is dominated by the full nonlinearity of the 
driving force but remains still regular. Within this do¬ 
main further bifurcations occur, however, without affect¬ 
ing the dc-current through the JJ which stays basically 
at its value somewhat above Ai. Phase-space orbits in¬ 
creasingly loose their harmonic-like ellipsoidal structure 
and become potato-shaped with deformations in form of 
‘dips’, thus reflecting the fact that the effective time- 
dependent potential I 4 ff((/?) ( 6 ) turns from mono-stable 
to multi-stable. 


(iii) With even further increasing A > A 2 (also Sec. IIIB), 
the system displays strong sensitivity to initial conditions 
(multiple steady-state orbits) and eventually irregular 
phase space patterns. When this domain is approached, 
the dc-current through the JJ grows substantially in par¬ 
allel with larger current variations depending on the driv¬ 
ing amplitude. 

(iv) For extremely strong driving A > A 3 ^ 1 and finite 
friction (Sec. IIIC), the dynamics turns into a regular 
motion again with a large time scale separation between 
the motion happening within one of the local wells in 
I 4 ff(t)<i 3 ) and global dynamics exploring I 4 ff over a wide 
range of tp. 


A. Weak driving regime 

As argued above, phase-space orbits in the weak driv¬ 
ing regime remain basically ellipses, while the amplitude 
becomes nonlinear until a bifurcation occurs. Thus, the 
steady-state orbit (7) is dominated by Fourier coefficients 
T’Oj <P±i- The corresponding equations are easily obtained 
by projecting Eq. ( 8 ) on the respective Fourier modes as 


ipo + AJi(</3i)cos((/3o - ^'i) = 0 


and with Ai = — I 

(pi[cos( 6 »i)Ai -b 7 nsin( 6 'i)] - A[Jo(v5i) sin((^o) - J2(7 ’i) sin((/?o - 26»i)] = 0 
ipi[- sin(0i)Ai -I- 7flcos(0i)] - A[Jo(</3i) cos{ipo) + J 2 (<Pi) cos((/?o - 20i)] = 0. 


( 10 a) 

( 10 b) 

( 10 c) 


This set of equations for {(/Jq, T’l, can be solved ap¬ 
proximately for small by exploiting that for the 

Bessel functions one has Jk{x) ~ 0{x^) for |a;| <i; 
1. One first gains from (10a) in leading order tpo ~ 
—AJi(i^i) cos( 0 i) « —Ay)icos( 0 i) so that driving depen¬ 
dent terms in ( 10 b) are of higher order in the small pa¬ 
rameter A. This yields the phase of oscillations, 

tan(6li) = -^ (11) 

with 01 = 0 at resonance. Inserting this result into (10c) 
leads close to resonance and in leading order in ipQ to 

o [’^o(v^i) + hipi)] ■ (12) 

a/A| -t 0272 

The known result for a driven harmonic oscillator is re¬ 
gained for \ipi\ 1 while for somewhat larger driving 
nonlinearities in the Bessel functions tend to play a role. 
This type of orbit, with amplitude named (p{ henceforth, 
which continually evolves from the harmonic-oscillator 
type of solution, is the only stable orbit until at a crit¬ 
ical driving strength Ai a second solution appears with 
amplitude . 


In contrast to the harmonic-oscillator type orbit 
the new orbit ip{^ exists even in absence of dissipation, 
and it is in this limit that it can easily be found analyt¬ 
ically: putting 7 = 0 in (10b), (10c) taken at resonance 
Ai = 0, we assume 933 = 0 to find the phase 9[^ = 7r/2 
from (10a). Its amplitude then follows via (10c) from 

To(</5i) — J2{pi) = 2-t— Ji{}Pi) = 0 (13) 

dip I 

as Piji — 1.841 independent of the driving amplitude. 
For finite friction and away from resonance orbits are 
obtained numerically from (4). 

It turns out that in a range beyond the threshold 
Ai the type-II orbit is the only stable one, while the 
harmonic-oscillator type-I orbits become unstable at this 
bifurcation point. The threshold Ai is determined by 
the condition that the amplitudes of both solutions 
match, :/3 ((Ai) = , i.e., according to (12) A 1/7 = 

« 2.912. In phase space both types of 
solutions display ellipsoidal orbits. The transition from 
type-I to type-II solutions at the bifurcation point Ai has 
also been found in the classical limit of a quantum de¬ 
scription within rotating-wave approximation in^®’^^. 
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Figure 3. (color online) Left: Mean dc-current (7j) through 
the JJ in steady state and on resonance = 1 as a function 
of the scaled driving amplitude A/y for (from bottom to top) 
7=0.005 (blue), 0.01 (black, cf. Fig. 2), 0.02 (red). 

Right: Power transfer Pjj^ho from the JJ to the resonator 
during one period of the oscillations in the steady state at 
too = n-2'K,n 1, n G N. The power transfer is shown for 7 = 

0.01 and at four different driving strengths (also marked in 
the left panel): A/ 7 = 0.1 (A), 2 (B), 3 (C) and 3.5 (D). While 
the harmonic-oscillator type-I solutions gain energy from the 
drive nearly during the whole oscillation cycle (A and B), 
for type-II oscillations the energy gained during part of the 
cycle is partly flowing back to the drive in other parts (C and 
D). In consequence, increasing the driving strength beyond 
A 1/7 « 2.9 does not further increase oscillation amplitude 
and current. 


An intuitive understanding of the saturation of the 
oscillation amplitude, when the driving strength is in¬ 
creased beyond Ai and, indeed, of the nature of the type- 
II orbit is offered by the numerical results in Fig. 3. How 
the resonator acts back onto the charge transfer in the 
JJ below/above the threshold is monitored by the en¬ 
ergy transferred from the JJ to the resonator [cf. (5)], 
i.e., the power Pjj_>.ho (0 = sin[(/:(t) -|- Ht] in the 

right panel of Fig. 3. Sufhciently below the threshold 
A 1/7 « 2.9, energy is nearly unidirectionally injected into 
the resonator, i.e., the drive pushes energy into the oscil¬ 
lator at each time of the oscillation cycle. The phase shift 
of the type-II oscillations, however, results in an oscilla¬ 
tion which extracts energy from the drive during one part 
of the oscillation cycle, and pushes energy back during 
another part. Increasing the driving strength beyond Ai 
will thus increase energy in- and back-flow, but will not 
further increase the net gain over a full cycle; hence, the 
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(color online) Top: Dissipated power Pdiaa from 
the driven JJ-|-resonator system for 7 = 0 . 01 , = 1 and 

various driving strengths Ai < A < A 2 around the second 
bifurcation Aj. Bottom: Snap-shots of the effective time- 
dependent potential Ves{ip) ( 6 ) for various driving strengths 
(see top for the color code) beyond the first threshold A > Ai 
and at times = (2n -|- l) 7 r, n € N. 
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saturation of the oscillation amplitude at — 1.841. 
Experimentally, the transition from type-I to type-II or¬ 
bits is seen as a saturation in the dc-current through the 
JJ, i.e., (Ij) = A(sin[(/j(t)-|-Ht])o , where the time average 
(•)n is taken in steady state and over several oscillation 
periods, see left panel of Fig. 3. The dc-current result¬ 
ing from this averaging is (approximately) proportional 
{(fl) [analytically found from Eq. (12) and (13)]. Physi¬ 
cally, this proportionality origins from balancing dissipa¬ 
tion (proportional to the stored energy) and power input. 
Accordingly, the quadratic dependence (Ij) oc oc Ej 
in the regime of very weak driving passes over to a lin¬ 
ear dependence for somewhat stronger driving but still 
before the threshold. This changeover is reminiscent of 
the changeover from the perturbative domain of sequen¬ 
tial (Coulomb blockade) to the regime of coherent charge 
transfer (phase coherent regime). Note that the various 
numerically found current curves in Fig. 3 approximately 
just scale with 7 , but for stronger damping the bifur¬ 
cation threshold is shifted below the result A 1/7 « 2.9 
calculated for the 7—^0 limit above. 


B. Beyond the first threshold 

Beyond the first threshold A > Ai, type-II orbits dic¬ 
tate the dynamics until a second bifurcation occurs at 
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Figure 5. (color online) Various classes of steady-state orbits (outer panels) depending on the initial conditions (inner panels) 
for 7 = 0.01 and strong A = 1 (left) and very strong driving A = 2.455 (right). The color code for the initial conditions in phase 
space corresponds to a particular type of steady-state orbit in phase space. Red dots result from Poincare plots and indicate 
points of return for steady-state orbits after one period 27r/n. 


AJ « 0.8 with only a very weak dependence on the fric¬ 
tion strength in this underdamped regime. The emer¬ 
gence of this class of orbits is characterized by a substan¬ 
tial deformation of the ellipsoidal phase-space structure 
(see Fig. 2), which is related to a changeover of the ef¬ 
fective potential T4ff((/j) (6) from being essentially mono¬ 
stable to dominantly multi-stable (see Fig. 4). 

A particle moving (rather weakly damped) in the time- 
dependent effective potential will encounter a deep well 
during its passage through </? = 0 in one direction, while 
on the way back in the second part of its oscillation cycle 
(see lower panel of Fig. 4), the well is less pronounced 
(corresponding to the situation below AJ) or even turns 
into a barrier around </? = 0 (situation well above A)(). In¬ 
deed, in the power dissipated from the driven system into 
the reservoir Pdiss{t) = —(upper panel of Fig. 4), 
the maximal amplitude at fit = 2mT (n integer) shows 
a local maximum right around the bifurcation value A*, 
whereas the speed of the fictitious particle (and concur¬ 
rently Pdiss) at nt = (2n + l)7r is more and more reduced 
and even develops a local minimum when passing through 
If = 0. While this second bifurcation has almost no effect 
on the dc-current (cf. Fig. 2), its appearance is detectable 
in the discussed features of the dissipated power. 

With further increasing the driving A > AJ and initial 
conditions close to the phase space origin, further bifur¬ 
cations occur that we do not need to discuss in detail 
here. It is important to note though that each bifurca¬ 
tion is associated with a change in stability meaning that 
only the newly emerging orbits determine the dynamics 


beyond each bifurcation threshold. Independent of the 
appearance of new orbits, the dc-current /j,dc stays ba¬ 
sically constant over a wide range of driving amplitudes 
Ai < A < A 2 « 1.6 (see Fig. 2). Similar to A* the numer¬ 
ical value of A 2 depends only very weakly on the friction 
strength in the underdamped regime. 

For the driving strengths considered so far, it has been 
sufficient to consider initial conditions {(/?(0), (/5(0)} close 
to the phase space origin only. Now, that steady-state 
orbits tend to explore larger domains in phase space, one 
may wonder about the impact of initial conditions lo¬ 
cated in these regions. This is illustrated in Fig. 5 (left), 
where classes of steady-state orbits are studied depending 
on their initial conditions for moderately strong driving 
A = 1, i.e., Ai < A < A 2 . 

Indeed, for different sets of initial conditions 
{y)(0), </5(0)}, one now asymptotically finds three different 
classes of steady-state orbits in phase space. This dynam¬ 
ical multi-stability reflects the multi-stability of the effec¬ 
tive potential I4ff(</?) due to a strong nonlinearity. How¬ 
ever, in this regime of driving-strengths, Ai < A < A 2 , 
the existence of multiple types of steady-state orbits will 
usually not be of experimental relevance. This is due to 
the fact, that multi-stability only occurs if the fictitious 
particle leaves local wells in I4ff which requires at least an 
energy of 2A . This could, in principle, be done by prepar¬ 
ing the system initially (a difficult task though) such that 
^'Ho(</2(0), </5(0)) > 2A: See the dark blue area in the 
left panel of Fig. 5, indicating initial conditions around 
the phase space center which lead to one and the same 
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steady-state orbit. Starting from an equilibrated circuit, 
however, domains of initial conditions leading to addi¬ 
tional steady-state orbits, do play a role only at highly 
elevated temperatures, ksT/muiQ > 2A. For the current 
experimental situation, this regime is not relevant. 

For A > A 2 « 1.6, one enters again a qualitatively 
new regime. It is characterized by a sharp rise of the dc- 
current through the JJ by almost an order of magnitude 
(see Fig. 2). Now, even for initial conditions close to the 
phase space origin, multiple types of stable steady-state 
orbits exist, both exploring large domains in phase space, 
see Fig. 5 (right). Some of these orbits are covered only 
after multiples of the fundamental period 27r/n, in con¬ 
trast to the regime A < A 2 [cf. single (multiple or single) 
Poincare-plot points in the orbits of the left (right) panel 
of Fig. 5]. Hence, the sensitivity with respect to initial 
conditions grows substantially, thus marking the onset of 
irregular, chaotic-like behavior for sufficiently large driv¬ 
ing amplitudes A > A 2 . A detailed analysis of properties 
and characteristics of possible chaotic dynamics in this 
domain is beyond the scope of this paper and will be 
presented elsewhere. 


C. Prom multi-well to elevator dynamics 

Interestingly, in the regime of extremely strong driv¬ 
ing, A > A 3 ^ 1, (and for finite friction) regular dynam- 





Figure 6 . (color online) (a) Phase space portraits of steady- 
state orbits with Poincare points (red) in the regime of ex¬ 
tremely strong driving with A = 30 > A 3 and 7 = 0.01 (top 
left), 0.025 (top right), 2.5 (bottom right), 5 (bottom left). 
Sketched in (b) is the particle dynamics in the effective poten¬ 
tial I4ff(v5), see Eq. ( 6 ), corresponding to phase space orbits in 
the right panels of (a). For moderate damping the particle is 
trapped in one of the local wells, ‘elevated’ upwards in Vefi(<p) 
[the motion towards large negative (/9-values with small ip in 
(a)] until its escape and consequent retrapping. 


ics dominates again (see Fig. 6 ). At this driving strength 
the periodic part of t 4 ff is pronounced enough to create a 
multitude of local minima in the superimposed quadratic 
potential, see Fig. 6 (b). We mention in passing that this 
multi-well pattern of I 4 ff has some analogy to the po¬ 
tential profile of superconducting quantum interference 
devices (SQUIDs). In contrast to SQUIDs, however, the 
potential Vee{^), cf. Eq. ( 6 ), combines a static quadratic 
with a time-dependent sinusoidal potential. 

The dynamics is then easily understood; most simply 
in the completely underdamped and the strongly over¬ 
damped cases. In the completely underdamped case [up¬ 
per left of Fig. 6 (a)] the particle essentially undergoes a 
simple oscillation in the quadratic potential over a wide ip 
region, running with high energy over the potential wells, 
which then causes slight wiggles in the phase-space orbit. 

For somewhat stronger friction [upper right of 
Fig. 6 (a)], a fictitious particle runs periodically through a 
cycle of localized and de-localized motion: Starting some¬ 
where in the low energy sector, it gets trapped quickly 
in one of the local minima of the potential close to the 
global minimum of Vesiip) ( 6 ). It is then transferred up 
in energy by the driving term oc cos ((/9 -I- Ht) while being 
trapped in this local well. During this process the poten¬ 
tial barrier of the respective local well shrinks until the 
particle can escape to run towards the global minimum 
while loosing substantial energy so that it gets trapped 
close to the global energy minimum again. This type of 
’’elevator” dynamics leads to increasingly simpler phase 
space patterns towards the overdamped regime. Phys¬ 
ically, during the trapping period jt^j grows almost lin¬ 
early with \ip\ ^ 1 /H, whereas the motion towards the 
global minimum is associated with an almost instanta¬ 
neous drop in amplitude accompanied by a large increase 
in momentum. As a consequence, based on the second 
Josephson relation V{t) cx ip{t), one expects to observe 
strong voltages pulses with frequencies much lower than 
the driving frequency D. 


IV. HIGHER ORDER RESONANCES 

So far we have considered driving the system at (or 
close to) the eigenfrequency, iu;o = 1 , of the resonator. 
At sufficiently strong driving, the system’s response then 
contains higher harmonics. It is, thus, also interesting 
to drive at frequencies D 9 ^ 1 , where higher or subhar¬ 
monics of the drive can become resonant with the eigen¬ 
frequency. Physically, the drive is detuned, of course, 
simply by changing the applied dc-voltage bias. 

Figure 7 shows the steady-state power spectrum, 
\(p(uj)\'^, i.e., the response of the system at frequency 
uj if driven at O. Such spectra have been recently in¬ 
vestigated experimentally in Ref. [29]. Resonances [i.e., 
stable steady-state solutions of (4) or ( 8 ) with large am¬ 
plitudes] are found for driving frequencies = n for in¬ 
teger n 7 ^ 0, where the system responds atw = 1,2,3,..., 
and for driving frequencies Hi = 1 /n with response fre- 
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Figure 7. (color online) a) Power spectrum, ln(|(^(a;)p), of 
steady-state orbits ip{t) for a range of driving frequencies Q 
in the underdamped regime 7 = 0.1 and at moderate driving 
amplitude A = 0.2. For the sake of presentation the spectrum 
is numerically broadened by taking the Fourier transform of a 
finite time signal, b) Dominant components of (p{t) [see 
(7)1, corresponding to lines uj = kQ in a). Panel c) shows a cut 
of the power spectrum at the (shifted) subharmonic resonance 
Q = 0.365 « fli at a larger driving amplitude A = 0.7, where 
the system dominantly responds with lo = 3D « 1 . 


quencies uj = 1/n, 2/n, 3/n,.... The situation, where the 
system is driven with D and responds with w < D is 
called down-conversion, the situation, where it responds 
with w > D up-conversion. Apparently, for a driving fre¬ 
quency D 2 the system dominantly responds with ui = 1 
while contributions of higher harmonics are weak. In con¬ 
trast, driving with Di one observes a response in which 
dominantly frequencies with w = 1/2 but also the first 
few higher harmonics ui = 1, 3/2,... are present. While 
the power spectrum, |(^(a;)p, shown here gives a direct 
intuitive link to the system dynamics, comparable in¬ 
formation can also be extracted from the directly mea¬ 
surable spectrum of light emission. Before we turn to 
details of its resonant features below, let us address the 
consequences of varying the dc-bias (and thus the driv¬ 
ing frequency) for the dc-current through the JJ. As de¬ 
picted in Fig. 8, current peaks appear indeed at driving 
frequencies and D « Du with n > 1 integer with a 
characteristic shift towards D > Di in the subharmonic 
domain. These features may remind of Shapiro steps^ 
known from JJs subject to dc- and ac-voltages of the 
form V{t) = Vdc + "Caccos(Dact), an interpretation which 
for the present situation is somewhat confusing though: 
The conventional argument for the existence of Shapiro 
steps is that resonances appear whenever the energy gap 
induced by the dc-voltage 2evdc matches multiples of ac- 
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Figure 8 . (color online) DC-current through the JJ vs. driv¬ 
ing frequency D and driving amplitude A for 7 = 0.1. For 
weak driving the system is nearly linear and dominated by 
the fundamental resonance, D = cuq = 1 , while for stronger 
driving the subharmonic resonances atD = Di_=wo/n and 
the multi-photon resonances at D = D„ = nujo become ap¬ 
parent. Note, the sharp onset of two-photon processes at 
Ac ~ 7 D = 0.2 described in Subsec. IV A, while subharmonic 
resonances increase smoothly and shift with increasing A as 
discussed in Subsec. IV B (cf. also Fig. 9). 


photon quanta hflac, ke., 2evdc/^ = nVlac, n > 1 being 
integer. For the present situation, however, one could 
argue in two ways (we temporarily return to physical di¬ 
mensions) : (i) The dc-voltage which determines the driv¬ 
ing frequency D = 2eV/Ti must be multiples of the res¬ 
onator frequency ojq implying D = nwo, n > 1 integer, or 
(ii) the excitation of the resonator by one energy quantum 
Tiujq requires multiples of ac-photon quanta nfiQ., n > 1 
integer, i.e., D = oj^/n. Along these lines, one could inter- 
prete either the subharmonic Dj, or the higher harmonic 
resonances as Shapiro steps. The fundamental differ¬ 
ence to the conventional set-up to observe Shapiro steps 
is that there the dynamics of the Josephson phase (j) is 
fixed by the external voltage according to 2eV{t) = ^/?i, 
while here the Josephson phase appears as a dynamical 
degree of freedom fixed by the dynamics of the resonator 
phase [see Eq. (2)]. 

A. Integer multi-photon processes 

Now, we start with multi-photon resonances = n, 
where the system asymptotically responds mostly with 
frequency a; = 1 (down-conversion). In particular, we 
focus on the generic case D = 2 (see Fig. 7) which can be 
interpreted as a parametric resonance. 

The general argument for the appearance of this type 
of resonances can be directly read off from the equation 
for the Fourier mode amplitudes (8): For D = n we seek 
orbits with a time dependence dominated by cos(t) im¬ 
plying that Fourier coefficients i with kn = 1 dominate 

the expansion (7). Accordingly, in (8) F± « Ff{(pi) and 
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upon projecting onto orbits with one arrives at time- 
independent equations for the coefficients Lpi which in¬ 
clude Bessel functions and J„+i(vJi) [see (9)]. 

In the particular case of = 2, we may write for small 
amplitudes A sin[(p(t)-|-2t] « A[sin(2t)-|-(/?(t) cos(2t)], thus 
giving in the equation of motion (4) effectively rise to a 
parametrically driven harmonic oscillator scenario. For 
this system steady-state solutions only exist below a crit¬ 
ical, damping dependent value Ac = 27 , while above that 
threshold oscillation amplitudes will grow infinitely in 
time. Of course, here the nonlinearity of the full problem 
prevents this divergence to occur. Then, the simplest 
ansatz, if{t) ~ (^1 cos(t -I- ), taking only oscillations at 

the oscillator’s resonance frequency into account, yields 
for the amplitudes 

+J 3 {^pi)\/ipi=-f/X. (14) 

Solutions of this equation only exist above the 
parametric-resonance threshold, X > Xc — 27 , with a 
phase 01 « — 7 r /4 (cf. also [36]). The full solution with 
an additional contribution <C and 9i « tt/2 
shows similar threshold behavior with some quantitative 
corrections (cf. the threshold in Fig. 8 and Fig. 9). 


B. Subharmonic resonances 

In the driving-frequency range < 1 below the fun¬ 
damental resonance, the system will dominantly respond 
with a few higher harmonics w = ft, 2ft,, see Fig. 7. 
Resonances appear approximately at driving frequencies 
fti, where one of these harmonics matches the eigenfre- 
quency wq = 1. This resonant response results, in par¬ 
ticular, in a strong contribution to the dc-current close 
to fli, where a shift of the resonance with increasing 
driving strength can be observed, see Fig. 8 . 

As a specific example, we consider the shift of the 
resonance close to 17 = 1/2 for 7 —0. Taking then 
an ansatz for the steady-state orbits (7) of the form 
ip{t) ~ ipi cos{ftt+ 6 i)+(p 2 cos{ 2 ftt+ 62 ) with 17 « IIi and 
considering ( 8 ) we find, that the Fourier coefficients for 
this type of steady-state orbits have to be gained from 
nonlinear equations including products of Bessel func¬ 
tions Assuming for suffi¬ 

ciently strong driving a dominant response with 217 « 1 
and thus (pi 1 , ip 2 , 'Vfe take into account only the low¬ 
est order of the Bessel functions for ipi and find the res¬ 
onance condition A 2 = —Xcos 0 i(piJi(ip 2 )/p >2 > 0 from 
the 2r7-projection of Eq. ( 8 ). That this shift is indeed 
positive, as seen in Fig. 8 , can be analytically confirmed 
for this limit by analyzing the 17-projection. Increasing 
the driving strength, the shift that depends on the mixing 
between the driving frequency and the first higher har¬ 
monics grows further. Remarkably, at least for moderate 
driving similar shifts do not occur for the r7„-resonances 
(including the fundamental one), as apparent from the 
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Figure 9. (color online) Time-averaged steady-state energy in 
the resonator for 7 = 0.1 and at different driving strengths: 
A =0.1 (black), 0.15 (blue), 0.2 (green), 0.5 (red). For the two- 
photon resonance, ft = 2 , there is a rather sharp threshold 
at A > Ac = 27 , while subharmonic resonances, ft ~ 1 /n 
gradually increase and shift with increased driving strength. 


numerical data in Fig. 8 and analytically following Eq. 
(10) and Sec. IV A. 

Note, that the multi-photon processes and subhar¬ 
monic resonances discussed here and shown in Figs. 7, 
8 all occur for comparatively weak driving; for strong 
driving similarly rich behavior for driving at sub- and 
higher-harmonic resonance frequencies may be expected, 
as found in Sec. Ill B, III C for the fundamental reso¬ 
nance. 


C. Energy transfer 

As already discussed for the fundamental resonance in 
Sec. Ill A, the energy transfer from the driving source 
to the resonator is another tool to reveal details of the 
nonlinear dynamics. Experimentally, it is accessible as 
(mean) photon emission from the resonator. Here, we 
depict resonance curves of the time averaged resonator 
energy Aho = -b :/ 3 ^/ 2 )o, see Eig. 9. As expected, 

the resonances discussed in the previous sections appear 
in form of pronounced peaks at frequencies f7„ and 17 1 . 
In the subharmonic regime the shift in the resonances 
towards 17 > f7i for larger A is seen as well. We note also 
the threshold A > Ac for the occurrence of the parametric 
resonance at 172 - 


V. THERMAL NOISE 

In actual experimental realizations, noise stemming 
from various sources is always present and may sensi¬ 
tively influence the dynamics of the JJ-bresonator device. 
In the Langevin equation (4) we restricted ourselves to 








10 



A 


Figure 10. (color online) Mean steady-state oscillation ampli¬ 
tude, i/Pmax = max[((p(t —^ oo))g], vs. driving strength aver¬ 
aged over 10000 realizations of thermal noise at temperatures 
ksT—O (black), 0.01 (green-dashed), 0.1 (blue), 0.25 (red) for 
friction parameter 7 = 0.01. 


thermal noise related to the finite photon lifetime in the 
resonator via the fluctuation-dissipation theorem. An¬ 
other major source of noise may be local voltage fluc¬ 
tuations at the JJ that may induce charge localization 
(Coulomb blockade). However, for the present situation 
the impact of these fluctuations is of minor relevance. Al¬ 
though assuming a purely classical regime implies, that 
temperatures are high as compared to quantum fluctua¬ 
tions, the circuit is nonetheless operated at temperatures 
sufficiently low compared to other energy scales of the 
system, i.e., in dimensional units k^T < mujQ,Ej. 

The simplest way to include finite temperature effects 
is to perform a thermal averaging over initial conditions, 
i.e., we start initially with a thermal distribution of phase 
space variables but then follow a deterministic time evo¬ 
lution according to (4) for ^ = 0. While somewhat in¬ 
consistent, this scenario accounts for the fact that precise 
initial conditions are experimentally not feasible. From 
a purely theoretical perspective, it allows to analyze the 
sensitivity of steady-state orbits onto initial conditions. 

The full description of thermal noise works with the 
Langevin equation (4) and seeks for phase space orbits 
averaged over many noise realizations {ip(t))^, {(p(t))^. In 
this latter situation, asymptotically thermal fluctuations 
may induce transitions between various steady-state or¬ 
bits even when their respective initial conditions are well 
separated in phase space. 

The main effect of an initial thermal distribution is 
apparent in the regime A > A 2 (see Fig. 2), where bare 
steady states tend to depend sensitively on the initial 
conditions. Accordingly, when averaged over an initial 
thermal distribution, phase space structures in steady 
state are washed out. This in turn gives rise to relatively 
large current fluctuations {{Idc — (ddc)o)^)o, where (• • • )o 



Figure 11. (color online) Mean dissipated power at temper¬ 
ature ksT = 0 (top two panels) and ksT = 0.1 (bottom) 
with 7 = 0.1 and = 2 for various driving strengths A=0.15 
(black), 0.2 (blue), 0.5 (green). The threshold for the para¬ 
metric amplihcation in absence of noise is Ac ~ 27 = 0.2. 


denotes the average over initial conditions according to a 
thermal distribution. In fact, the size of current fluctua¬ 
tions directly indicates ranges in parameter space, where 
the underlying asymptotic dynamics displays either bi¬ 
furcations or irregular behavior, cf. Fig. (2). 

Thermal noise according to the full Langevin dynam¬ 
ics has a similar impact. Here, we focus on two do¬ 
mains, namely, the domain around the first bifurcation 
A « Ai for fl = 1 (see Sec. HI A) and the domain around 
the parametric resonance A « Ac « 27 for H = 2 (see 
Sec. IV A). 

For the first case, in Fig. 10 the mean steady-state am¬ 
plitude max [{ip{t —>■ 00 ))^] is shown. Even in the weak 
driving regime the linear response of the system gets in¬ 
fluenced by temperature as thermal fluctuations become 
larger and increasingly explore the nonlinearity of the ef¬ 
fective potential. More substantial deviations occur for 
A —>■ Ai, where the bifurcation is increasingly smeared 
out and the overall oscillation amplitude decreases at el¬ 
evated temperatures. 

In the second case, the parametric-resonance thresh¬ 
old, the situation is a bit more intricate. Here, thermal 
noise may mix higher harmonics in the Fourier expansion 
(7) with the consequence that down-conversion already 
occurs prior to the threshold Ac, see Fig. 11. Accordingly, 
already for driving strengths below the zero-temperature 
threshold, thermal noise excites oscillations with the res¬ 
onant frequency a; = 172/2 = 1 with a drastically in¬ 
creased amplitude compared to the response at the driv¬ 
ing frequency w = 172 in the absence of noise (e.g., an 
increase by about two orders of magnitude for the black 
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lines in Fig. 11). Above-threshold oscillations are some¬ 
what reduced by temperature. Nonetheless, a transition 
remains clearly visible, for example, in the mean dissi¬ 
pated power (Pdiss(t)){- Furthermore, an offset emerges 
such that (F’diss(t)){ > 0 for all times, which can be re¬ 
lated to the thermal energy continuously injected into the 
system. 

VI. CONCLUSION AND DISCUSSION 

In this paper we analyzed the classical dynamics of a 
circuit, where a single relevant resonator mode interacts 
with a dc-voltage biased JJ. This problem can be mapped 
onto the dissipative dynamics of a fictitious particle mov¬ 
ing in a nonlinear, time dependent potential, where in 
contrast to conventional settings the nonlinearity appears 
as part of the driving, while the static part of the po¬ 
tential is purely harmonic. In the regime of moderate to 
large Q-factors (underdamped regime) and weak thermal 
noise, steady-state orbits and corresponding observables 
are determined by basically only two parameters, namely, 
the dimensionless driving strength (Josephson coupling) 
A = Ej /towq and the dimensionless driving frequency fl 
(in units of wq)- 

At the fundamental resonance fl = tuj jujQ = 1 this sys¬ 
tem displays a changeover from a linear response regime 
for weak driving towards a strongly nonlinear behav¬ 
ior for strong driving. The various dynamical domains 
leave their signatures in the dc-current flowing through 
the JJ and in the microwave power emitted from the 
resonator and are thus directly accessible experimen¬ 
tally. Resonances are also found when driving with either 
higher harmonics = n, n integer) or sub-harmonics 
(n « l/u, n integer), while the system responds with 
the fundamental frequency thus corresponding to pro¬ 
cesses of down- and up conversion. These features can 
also be detected by either monitoring the dc-Josephson 
current or the radiated microwaves. Due to its high de¬ 
gree of tunability the resonator-l-JJ circuit thus allows 
to study the full wealth of classical nonlinear dynamics 
in one-dimensional driven, dissipative systems. The im¬ 
pact of weak thermal noise is most prominent close to 
bifurcations of steady-state orbits. 

At this point let us discuss a typical set of parame¬ 
ters for circuit designs that allows to access the physics 
discussed above: The classical regime with weak ther¬ 


mal noise requires that towq ^ k^T ^ hujQ. For a res¬ 
onance frequency of ojq ^ 5 GHz, this can be realized 
with an LC-circuit with (7^5 pF operated at temper¬ 
atures T ^ 150 mK. One then has mujQ ~ 0.08 meV, 
k^T ~ 0.015 meV, and Tiujq ^ 0.003 meV. In the phase 
regime of the JJ, one further needs Ej ^ Eq which 
applies for a typical Ej ~ 1 meV. An external mag¬ 
netic flux allows to tune this maximal Josephson cou¬ 
pling down to a factor of about 100 , i.e., within the 
range Ej ^ 0.01... 1 meV. Present resonator designs 
have typical photon lifetimes over a wide range of Q- 
factors Q « 10 ... 10 '* which coincides with the (strongly) 
underdamped regime. 

Theoretically, at the fundamental resonance new dy¬ 
namical domains are associated with driving parame¬ 
ters Ai,A*,A 2 , and A 3 as depicted in Fig. 2 which, given 
the above parameters, translates into the following cou¬ 
pling energies: Ej^i « 0.02 meV, Ej^^ « 0.05 meV, 
Ej ^2 ~ 0.1 meV, Ej^ 3 {X « 20) « 1.3meV. Apparently, 
the first three coupling energies are easily accessible ex¬ 
perimentally with only Aya lying at the edge of the range. 
The challenge here is that to cover the full range of driv¬ 
ing amplitudes within one set-up from the linear response 
regime Ej <C i?yi to Ej^s, requires to vary Ej by a factor 
of somewhat more than 100 . 

Note that the superconducting gap 2A defines an up¬ 
per limit on the applied voltage 2eU < 2A which for A1 
junctions corresponds to wj ^ 500 GHz, thus allowing 
also for driving at higher harmonics Following this 
discussion, we are confident that the classical, nonlinear 
dynamics analyzed in this work is indeed accessible in re¬ 
alistic circuits similar to those in reference [29] and opens 
the door to study the interplay of driving, dissipation and 
resonances under well controlled conditions. 
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